
// 用于得到一系列剂量值的放射性尘埃等值线的各点、以及覆盖面积
class SFSS{
    constructor(){
        this.W = 0 ;
        this.fraction = 1;
        this.windvelocity = 0;
        this.windangle = 0;
        this.Is = [];
        this.contours_data = [];
        this.contours_point = [];
        this.gzlng = 0; // 爆心地面投影点经度坐标
        this.gzlat = 0; // 爆心地面投影点纬度坐标
        this.mushroomsize = {
            h: 0,
            a: 0,
            b: 0,
            r_s: 0
        };

    }

    init(burstyield, fraction, windvelocity, windangle, doses, gzlng, gzlat){
        this.contours_data = [];
        this.contours_point = [];
        this.W = burstyield;
        if((this.W < 1 || this.W > 1e5) && debug){
            console.log("SFSS: W out of range 1 ~ 1e5, W:", this.W);
        }
        this.fraction = fraction;
        this.windvelocity = windvelocity;
        this.windangle = windangle;
        this.Is = doses;
        this.gzlng = gzlng; // 爆心地面投影点经度坐标
        this.gzlat = gzlat; // 爆心地面投影点纬度坐标

        this.getContourParam();
        this.getPoints();
    }

    // 返回用于确定stem、cloud的参数
    getParam(){
        let a_2_3 = unlog10(-0.509 + 0.076 * log10(this.W));
        let a_4 = unlog10(0.270 + 0.089 * log10(this.W));
        let a_5 = unlog10(-0.176 + 0.022 * log10(this.W));
        let a_5_i = unlog10(-0.054 + 0.095 * log10(this.W));
        let a_6 = unlog10(0.030 + 0.036 * log10(this.W));
        let a_7 = unlog10(0.043 +0.141 * log10(this.W));
        let a_8 = unlog10(0.185 + 0.151 * log10(this.W));
        let a_9;
        if(1<= this.W && this.W < 28)
            a_9 = unlog10(1.371 - 0.124 * log10(this.W));
        else if (28 <= this.W && this.W <= 1e5)
            a_9 = unlog10(0.980 + 0.146 * log10(this.W));
        

        let a_s = unlog10(2.880 + 0.348 * log10(this.W));
        let a = unlog10(3.389 + 0.431 * log10(this.W));
        let b = 1.40* 1e3 *Math.pow(this.W, 0.431);
        let h;
        if(1<= this.W && this.W < 28)
            h = unlog10(3.820 + 0.445 * log10(this.W));
        else if (28 <= this.W && this.W <= 1e5)
            h = unlog10(4.226 + 0.164 * log10(this.W));
        let  R_s = unlog10(2.319 + 0.333 * log10(this.W)); 
        let log_a_Rs = 1.070 + 0.098 * log10(this.W);
        let a_o = unlog10(log10(a) - (h * log_a_Rs)/(h - R_s));
        let k_a = 2.303 * (log_a_Rs/(h - R_s));
        let z_o;
        let z_s = 2.303 * (log10(a_s / a_o)) / k_a;
        if(1<= this.W && this.W < 9)
            z_o = h - b;
        else if (9 <= this.W && this.W <= 1e5)
            z_o = (1900 + (a_2_3 + 0.020) * z_s) / a_2_3;
        
        let X_1;
        let X_2;
        let X_3;
        let X_4;
        let X_5;
        let X_6;
        let X_7;
        let X_8;
        let X_9;
        X_2 = a_2_3 * z_s - a_s;
        X_3 = a_2_3 * z_s + a_s;
        X_4 = (a_4 * (a_4 * z_o - 1900)) / (a_4 + 0.020);

        if(1 <= this.W && this.W < 28){
            X_1 = - unlog10(3.308 + 0.496 * log10(this.W));
            X_5 = unlog10(3.644 + 0.467 * log10(this.W));
            X_6 = unlog10(3.850 + 0.481 * log10(this.W));
            X_7 = unlog10(3.862 + 0.586 * log10(this.W));
            X_8 = unlog10(4.005 + 0.596 * log10(this.W));
            X_9 = unlog10(5.190 + 0.319 * log10(this.W))
        }
        else if(28 <= this.W && this.W <= 1e5){
            X_1 = - unlog10(3.564 + 0.319 * log10(this.W));
            X_5 = unlog10(4.049 + 0.186 * log10(this.W));
            X_6 = unlog10(4.255 + 0.200 * log10(this.W));
            X_7 = unlog10(4.268 + 0.305 * log10(this.W));
            X_8 = unlog10(4.410 + 0.315 * log10(this.W));
            X_9 = unlog10(5.202 + 0.311 * log10(this.W));
        }
        let Y_s = unlog10(3.223 + 0.400 * log10(this.W));

        let logWs = [0, 1, 2, 3, 4, 5];
        let Y_8s = [6620, 12200, 48200, 167000, 342000, 650000];
        let Y_8_index = logWs.findIndex(nolessthan, log10(this.W));
        let Y_8 = interpolate(logWs[Y_8_index - 1], Y_8s[Y_8_index - 1], logWs[Y_8_index], Y_8s[Y_8_index], log10(this.W));

        let I_1;
        let I_2_3;
        let I_4;
        let I_5;
        let I_6;
        let I_7;
        let I_8;
        let I_9;
        I_1 = I_4 = 1; // 1 rads/h
        let k_1_2;
        if(1 <= this.W && this.W < 28)
            // 2.503 前应该有-，参考自NUKEMAP
            k_1_2 = unlog10(-2.503 - 0.404 * log10(this.W));
        else if(28 <= this.W && this.W <= 1e5)
            k_1_2 = unlog10(-2.600 - 0.337 * log10(this.W));
        I_2_3 = unlog10(k_1_2 * (X_2 - X_1) / 2.303);
        let as = [];
        let a_h;
        if(1 <= this.W && this.W < 28)
            a_h = unlog10(-0.431 - 0.014 * log10(this.W));
        else if(28 <= this.W && this.W <= 1e5)
            a_h = unlog10(-0.837 + 0.267 * log10(this.W));
        let a_b_2 = unlog10(0.486 + 0.262 * log10(this.W));
        as[5] = a_5;
        as[6] = a_6;
        as[7] = a_7;
        as[8] = a_8;
        as[9] = a_9;
        let fais = [];
        for(let i = 5; i <= 9; i ++){
            if(as[i] >= a_h){
                fais[i] = (as[i] + a_h + Math.pow(a_b_2 + Math.pow(as[i] + a_h, 2), 1/2))
                          /(as[i] - a_h + Math.pow(a_b_2 + Math.pow(as[i] - a_h, 2), 1/2));
            }
            else if(as[i] < a_h){
                // 参考自NUKEMAP，分母的乘号应该为+
                fais[i] = (as[i] + a_h + Math.pow(a_b_2 + Math.pow(as[i] + a_h, 2), 1/2))
                         / (a_2_3 + Math.pow(a_b_2 + a_2_3 * a_2_3, 1/2));
            }
        }
        let kA_5;
        let kA_6;
        let kA_7;
        let kA_8;
        let kA_9;
        if(as[5] >= a_h){
            if(1 <= this.W && this.W < 28)
                kA_5 = unlog10(-3.286 - 0.298 * log10(this.W));
            else if(28 <= this.W && this.W <= 1e5){
                kA_5 = unlog10(-2.889 - 0.572 * log10(this.W));
            }
        }
        else if(as[5] < a_h){
            kA_5 = unlog10(-3.185 - 0.406 * log10(this.W));
        }
        if(as[6] >= a_h){
            kA_6 = unlog10(-1.134 - 0.074 * log10(this.W));
        }
        else if(as[6] < a_h){
            kA_6 = unlog10(-1.225 - 0.022 * log10(this.W));
        }
        if(as[7] >= a_h){
            kA_7 = unlog10(-0.989 - 0.037 * log10(this.W));
        }
        else if(as[7] < a_h){
            kA_7 = unlog10(-1.079 - 0.020 * log10(this.W));
        }
        if(as[9] >= a_h){
            kA_9 = unlog10(-2.166 - 0.552 * log10(this.W));
        }
        else if(as[9] < a_h){
            kA_9 = unlog10(-2.166 - 0.552 * log10(this.W));
        }
        I_5 = 4.606 * a * kA_5 * log10(fais[5]);
        I_6 = 4.606 * a * kA_6 * log10(fais[6]);
        I_7 = 4.606 * a * kA_7 * log10(fais[7]);
        I_9 = 4.606 * a * kA_9 * log10(fais[9]);

        this.mushroomsize.h = h;
        this.mushroomsize.a = a;
        this.mushroomsize.b = b;
        this.mushroomsize.r_s = R_s;
        
        return {
            X_1: X_1,
            X_2: X_2,
            X_3: X_3,
            X_4: X_4,
            X_5: X_5,
            X_6: X_6,
            X_7: X_7,
            X_8: X_8,
            X_9: X_9,
            Y_s: Y_s,
            Y_8: Y_8,
            I_1: I_1,
            I_2_3: I_2_3,
            I_4: I_4,
            I_5: I_5,
            I_6: I_6,
            I_7: I_7,
            I_9: I_9,
            z_o: z_o
        }
    }

    // 计算得到确定stem、cloud的形状的参数
    getContourParam(){
        let pm = this.getParam();
        if(debug) console.log('SFSS,pm', pm);
        for(var i = 0; i < this.Is.length; i ++){
            /** 
             contour_data包括：
                rads_per_h;
                draw_stem; 
                stem_upwinddistance;
                stem_downwinddistance;
                stem_maxY;
                stem_area;
                draw_cloud;
                cloud_upwinddistance;
                cloud_downwinddistance;
                cloud_maxY;
                cloud_area;
            */
            let contour_data = {};
            let stem_Xi_1;
            let stem_Xi_2;
            let cloud_X8i;
            let cloud_Xi_1;
            let cloud_Xi_2;
            let I =  this.Is[i] * 1 / this.fraction; // 裂变份额修正。对于裂变份额小于1的情况，绘制剂量率为原剂量率放大1/fraction倍的等值线

            contour_data.rads_per_h = this.Is[i];

           
            contour_data.draw_stem = I < pm.I_2_3;
            
            contour_data.draw_cloud = I < pm.I_7;
            
            if(contour_data.draw_stem){
                stem_Xi_1 = interpolate(log10(pm.I_1), pm.X_1, log10(pm.I_2_3), pm.X_2, log10(I));
                stem_Xi_2 = interpolate(log10(pm.I_2_3), pm.X_3, log10(pm.I_4), pm.X_4, log10(I));
                contour_data.stem_maxY = interpolate(log10(pm.I_2_3), 0, 0, pm.Y_s, log10(I)) * feet2m;
                contour_data.stem_upwinddistance = (pm.X_2 - stem_Xi_1) * feet2m;
                contour_data.stem_downwinddistance = (stem_Xi_2 - pm.X_2) * feet2m;
                contour_data.stem_center = pm.X_2 * feet2m;
            }
            else{
                contour_data.stem_maxY = 0.;
                contour_data.stem_upwinddistance = 0.;
                contour_data.stem_downwinddistance = 0.;
                contour_data.stem_center = 0.;
            }

            if(contour_data.draw_cloud){
                if(I <= pm.I_6)
                    cloud_Xi_1 = interpolate(log10(pm.I_5), pm.X_5, log10(pm.I_6), pm.X_6, log10(I));
                else if(I > pm.I_6)
                    cloud_Xi_1 = interpolate(log10(pm.I_6), pm.X_6, log10(pm.I_7), pm.X_7, log10(I));
                cloud_Xi_2 = interpolate(log10(pm.I_7), pm.X_7, log10(pm.I_9), pm.X_9, log10(I));
                cloud_X8i = interpolate(log10(pm.I_7), pm.X_7, 0, pm.X_8, log10(I));
                // cloud_maxY 用以下方式，与Fallout Models and Radiological Countermeasure Evaluations p137 Fg.21对照发现区别很大
                // contour_data.cloud_maxY = cloud_X8i * pm.Y_8 / (pm.X_8 - pm.X_7) * feet2m;
                // 采取插值：
                contour_data.cloud_maxY = pm.Y_8 * log10(pm.I_7/ I) / log10(pm.I_7 - pm.I_9);
                contour_data.cloud_upwinddistance = (cloud_X8i - cloud_Xi_1) * feet2m;
                contour_data.cloud_downwinddistance = (cloud_Xi_2 - cloud_X8i) * feet2m;
                contour_data.cloud_center = cloud_X8i * feet2m;
            }
            else{
                contour_data.cloud_center = 0.;
                contour_data.cloud_maxY = 0.;
                contour_data.cloud_upwinddistance = 0.;
                contour_data.cloud_downwinddistance = 0.;
            }

            // 风速修正
            let F;
            if(this.windvelocity > 15)
                F = 1 + (this.windvelocity - 15) / 60;
            else if(this.windvelocity <= 15)
                F = 1 + (this.windvelocity -15) / 30;
            contour_data.stem_downwinddistance *= F;
            contour_data.cloud_downwinddistance *= F;

            // // 面积
            // contour_data.cloud_area = Math.PI * contour_data.cloud_upwinddistance * contour_data.cloud_maxY / 2 + Math.PI * contour_data.cloud_downwinddistance * contour_data.cloud_maxY / 2;
            // contour_data.stem_area = 

            this.contours_data.push(contour_data);
        }
    }

    // 返回用于绘制等值线的点
    getPoints(){
        for(let i = 0; i < this.Is.length; i ++){
            let data = this.contours_data[i];
            let c_points = [];

            // 得到各点
            if(data.draw_stem){
                // stem 与 cloud 之间没有重合
                if((data.stem_center + data.stem_downwinddistance <= data.cloud_center - data.cloud_upwinddistance) && data.draw_cloud){
                    this.ellipseArc(c_points, data.stem_center - data.stem_upwinddistance, 0, data.stem_center, data.stem_upwinddistance);
                    let p1 = [];
                    let p2 = [];
                    this.ellipseArc(p1, data.stem_center, data.stem_upwinddistance, data.stem_center + data.stem_upwinddistance, 0);
                    // 简便起见，stem的右边的椭圆部分不考虑左上的1/4圆弧
                    this.ellipseArc(p2, data.stem_center, data.stem_maxY, data.stem_center + data.stem_downwinddistance, 0);
                    this.delInner(p1, data.stem_center, 0, data.stem_downwinddistance, data.stem_maxY);
                    this.delInner(p2, data.stem_center, 0, data.stem_upwinddistance, data.stem_upwinddistance);
                    this.adds(c_points, p1);
                    this.adds(c_points, p2);
                    
                    this.ellipseArc(c_points, data.cloud_center - data.cloud_upwinddistance, 0, data.cloud_center, data.cloud_maxY);
                    this.ellipseArc(c_points, data.cloud_center, data.cloud_maxY, data.cloud_center + data.cloud_downwinddistance, 0);
                    
                    // 关于X轴对称
                    for(let i = c_points.length - 1; i > -1; i --)
                        c_points.push([c_points[i][0], - c_points[i][1]]);
                }
                // stem 与 cloud 之间有重合
                else if((data.stem_center + data.stem_downwinddistance > data.cloud_center - data.cloud_upwinddistance) && data.draw_cloud){
                    this.ellipseArc(c_points, data.stem_center - data.stem_upwinddistance, 0, data.stem_center, data.stem_upwinddistance);
                    let p1 = [];
                    let p2 = [];
                    let p3 = [];
                    this.ellipseArc(p1, data.stem_center, data.stem_upwinddistance, data.stem_center + data.stem_upwinddistance, 0);
                    // 简便起见，stem的右边的椭圆部分不考虑左上的1/4圆弧
                    this.ellipseArc(p2, data.stem_center, data.stem_maxY, data.stem_center + data.stem_downwinddistance, 0);
                    this.ellipseArc(p3, data.cloud_center - data.cloud_upwinddistance, 0, data.cloud_center, data.cloud_maxY);
                    this.delInner(p1, data.stem_center, 0, data.stem_downwinddistance, data.stem_maxY);
                    this.delInner(p2, data.stem_center, 0, data.stem_upwinddistance, data.stem_upwinddistance);
                    this.delInner(p2, data.cloud_center, 0, data.cloud_upwinddistance, data.cloud_maxY);
                    this.delInner(p3, data.stem_center, 0, data.stem_downwinddistance, data.stem_maxY);
                    this.adds(c_points, p1);
                    this.adds(c_points, p2);
                    this.adds(c_points, p3);
                    this.ellipseArc(c_points, data.cloud_center, data.cloud_maxY, data.cloud_center + data.cloud_downwinddistance, 0);

                    // 关于X轴对称
                    for(let i = c_points.length - 1; i > -1; i --)
                        c_points.push([c_points[i][0], - c_points[i][1]]);
                }
                else if(!data.draw_cloud){
                    this.ellipseArc(c_points, data.stem_center - data.stem_upwinddistance, 0, data.stem_center, data.stem_upwinddistance);
                    let p1 = [];
                    let p2 = [];
                    this.ellipseArc(p1, data.stem_center, data.stem_upwinddistance, data.stem_center + data.stem_upwinddistance, 0);
                    // 简便起见，stem的右边的椭圆部分不考虑左上的1/4圆弧
                    this.ellipseArc(p2, data.stem_center, data.stem_maxY, data.stem_center + data.stem_downwinddistance, 0);
                    this.delInner(p1, data.stem_center, 0, data.stem_downwinddistance, data.stem_maxY);
                    this.delInner(p2, data.stem_center, 0, data.stem_upwinddistance, data.stem_upwinddistance);
                    this.adds(c_points, p1);
                    this.adds(c_points, p2);

                    // 关于X轴对称
                    for(let i = c_points.length - 1; i > -1; i --)
                        c_points.push([c_points[i][0], - c_points[i][1]]);
                }
            }
            else if(data.draw_cloud){
                this.ellipseArc(c_points, data.cloud_center - data.cloud_upwinddistance, 0, data.cloud_center, data.cloud_maxY);
                this.ellipseArc(c_points, data.cloud_center, data.cloud_maxY, data.cloud_center + data.cloud_downwinddistance, 0);
                // 关于X轴对称
                for(let i = c_points.length - 1; i > -1; i --)
                    c_points.push([c_points[i][0], - c_points[i][1]]);
            }
            else{
                if(debug)
                    console.log('SFSS.getPoints: 不能绘制出该等值线, doserate:', data.rads_per_h);
            }

            // 旋转各点
            this.rotate(c_points, this.windangle / 180 * Math.PI);

            // console.log('points m: p', c_points);
            // 转换为经纬度坐标, BMApGL Point 类型
            this.transToLngLat(c_points);

            this.contours_point.push(c_points);
        }
    }

    // 1/4椭圆弧各点。(x1, y1), (x2, y2)按逆时针方向取
    ellipseArc(p, x1, y1, x2, y2){
        let xc;
        let yc;
        let startAngel;
        let stopAngel;
        let a;
        let b;
        // 确定中心坐标, 起始角度和终止角度
        // 根据(x1, y1), (x2, y2)之间的关系分四种情况确定
        if(x1 < x2){
            if(y1 < y2){
                xc = x2;
                yc = y1;
                startAngel = Math.PI;
                stopAngel = Math.PI / 2;
            }
            else if(y1 > y2){
                xc = x1;
                yc = y2;
                startAngel = Math.PI / 2;
                stopAngel = 0;
            }
        }
        else if(x1 > x2){
            if(y1 > y2){
                xc = x2;
                yc = y1;
                startAngel = 0;
                stopAngel = - Math.PI / 2;
            }
            else if(y1 < y2){
                xc = x1;
                yc = y2;
                startAngel = - Math.PI / 2;
                stopAngel = - Math.PI;
            }    
        }
        
        // 半长轴a, 半短轴b
        a = Math.abs(x1 - x2);
        b = Math.abs(y1 - y2);

        // 得到点
        for(let angle = startAngel; angle >= stopAngel; angle += (stopAngel - startAngel)/100.){
            let x = xc + a * Math.cos(angle);
            let y = yc + b * Math.sin(angle);
            p.push([x, y]);
        }
    }

    // 判断点是否在椭圆内
    inEllipse(x, y, xc, yc, a, b){
        let isIn = (Math.pow(x - xc, 2) / Math.pow(a, 2) + Math.pow(y - yc, 2) / Math.pow(b, 2)) < 1;
        return isIn;
    }

    // 去除在内部的点
    delInner(p, xc, yc, a, b){
        let startIndex;
        let endIndex;
        for(let i = 0; i < p.length; i ++){
            if(this.inEllipse(p[i][0], p[i][1], xc, yc, a, b)){
                if(startIndex == undefined)
                    startIndex = i;
                endIndex = i;
            }     
        }
        let l = endIndex - startIndex + 1;
        if(startIndex == 0){
            for(let i = 0; i < l; i ++)
                p.shift();
        }
        else if(endIndex == (p.length - 1)){
            for(let i = 0; i < l; i ++)
                p.pop();
        }
    }

    // 将数组p2添加到数组p1末尾
    adds(p1, p2){
        for(let i = 0; i < p2.length; i ++){
            p1.push(p2[i]);
        }
    }

    // 旋转坐标
    rotate(p, angle){
        for(let i = 0; i < p.length; i ++){
            let x = p[i][0];
            let y = p[i][1];
            p[i][0] = x * Math.cos(angle) - y * Math.sin(angle);
            p[i][1] = x * Math.sin(angle) + y * Math.cos(angle);
        }
    }

    // 转换为经纬度坐标
    transToLngLat(p){
        // 假设地球为球体，将等值线各点坐标转换为经纬度坐标
        // 百度地图坐标先经度，再纬度。故：x -> lng, y -> lat, (x, y) -> (lng, lat)
        // 存在问题：利用百度地图纬度坐标得到所处纬度圆的半径是不准确的
        let R = 6371000;
        for(let i = 0; i < p.length; i++){
            let lng = this.gzlng + (p[i][0] / (R * Math.cos(this.gzlat * Math.PI / 180))) * 180 / Math.PI;
            let lat = this.gzlat + (p[i][1] / R) * 180 / Math.PI;
            // p[i][0] = lng;
            // p[i][1] = lat;
            p[i] = new BMap.Point(lng, lat);
        }
    }
}

if (typeof register == 'function') { register("fallout.js"); }